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The equilibrium behavior of vortices in the classical two-dimensional (2D) XY model with un- 
corrected random phase shifts is investigated. The model describes Josephson- Junction arrays with 
positional disorder, and has ramifications in a number of other bond-disordered 2D systems. The 
vortex Hamiltonian is that of a Coulomb gas in a background of quenched random dipoles, which 
is capable of forming either a dielectric insulator or a plasma. We confirm a recent suggestion by 
Nattermann, Scheidl, Korshunov, and Li [J. Phys. I (France) 5, 565 (1995)], and by Cha and Fertig 
[Phys. Rev. Lett. 74, 4867 (1995)] that, when the variance a of random phase shifts is smaller 
than a critical value er c , the system is in a phase with quasi-long-range order at low temperatures, 
without a reentrance transition. This conclusion is reached through a nearly exact calculation of the 
single- vortex free energy, and a Kosterlitz-type renormalization group analysis of screening and ran- 
dom polarization effects from vortex-antivortex pairs. The critical strength of disorder a c is found 
not to be universal, but generally lies in the range < a c < tv/8. Argument is presented to suggest 
that the system at a > er c does not possess long-range glassy order at any finite temperature. In 
the ordered phase, vortex pairs undergo a series of spatial and angular localization processes as the 
temperature is lowered. This behavior, which is common to many glass-forming systems, can be 
quantified through approximate mappings to the random energy model and to the directed polymer 
on the Cayley tree. Various critical properties at the order-disorder transition are calculated. 
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I. INTRODUCTION 



(^■)=0, {A%) = a, 
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The Kosterlitz-Thouless-Berezinskii (KTB) transition 
plays an important role in the theory of ordering 
in two-dimensional (2D) systems which has a continuous 
symmetry specified by a phase. Examples include planar 
magnets, 2D solids, Josephson- Junction arrays, and su- 
perfluid and superconductor films, etc. These systems 
have an ordered phase at low temperatures, characterized 
by power-law decay of correlations with distance. The 
(quasi)long-range order is destroyed through unbinding 
of vortex-antivortex pairs, which takes place at the KTB 
transition. 

A question of both theoretical and practical interest is 
whether and how quenched disorder alters the above pic- 
ture. In this paper we shall focus on the case of random 
frustration, where disorder introduces random, uncorre- 
cted phase shifts but do not pin the phase angles them- 
selves. More precisely, we shall consider an XY model 
with the following Hamiltonian Q , 



#({&}) = -J 5>os( 

(ij) 



(1.1) 



where the sum runs over all nearest neighbor pairs on 
a square lattice. The quenched random variables Ay, 
which give a random bias to the preferred advancing an- 
gle over each bond, are assumed to be uncorrelated from 
bond to bond, and each is gaussian distributed with the 
mean and variance given by 



respectively. It has been suggested that model ([O]) pro- 
vides a good description of the Josephson- Junction arrays 
in a transverse magnetic field [p[-p!o[. In this case, 4n is 
identified with the phase of the superconducting order 
parameter of grain i, and Ay = (27r/<I>o) A cxt • dl, 
where A oxt is the vector potential of the external mag- 
netic field and <E>o = h c/2 e is the superconducting flux 
quantum. The case (L2) corresponds to a situation 
where the average magnetic flux over each elementary 
plaquette of the grain network is an integer multiple of 
$0, but random displacement of superconducting grains 
from a perfect lattice structure yields quenched random 
phase shifts . 



On the theoretical side, model (1.1) and its variants 
have been studied extensively in the past J5|,[b2| [L9| . Re- 
sult of previous studies can be summarized as follows, 
(i) The spin-wave fluctuations have the same excitation 
spectrum as in the pure case. Disorder introduces dis- 
tortion in the ground state away from a perfect ferro- 
magnetic alignment. The combined effect of thermal and 
disorder fluctuations leads to an algebraic decay of the 
two-point phase-phase correlation function, 

C s „(r y ) = (exp[z(0 SWj4 - 0swj)]) ~ ry JJsw , (1.3) 

where m is the distance between site i and j, and 



(1.4) 
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is the correlation length exponent at temperature T, due 
to spin-waves only, (ii) Vortices, which are topological 
point defects in the 0- field, interact with each other and 
with the quenched disorder through a Coulomb potential. 
The interaction between two vortices is of the charge- 
charge type, where the charge of each vortex is given by 
its vorticity. The interaction between a vortex and a par- 
ticular disordered bond is of the charge-dipole type, with 
the strength of the dipole given by the phase shift Ay- 
over the bond. The equilibrium statistics of vortices is 
decoupled from that of spin waves. For a long time, the 
phase diagram of the model was thought to be of the 
type illustrated in Fig. 1(a) pj|jL^]. For a < a c ~ n/8, a 
phase with bound vortex-antivortex pairs, and hence al- 
gebraic decay of phase correlations, still exists, but only 
in a temperature window TL(cr) < T < T + (a). Below 
T_(cr), a "re-entrant" disordered phase was predicted. 
The two transition temperatures coincide at a critical 
strength of the disorder a c , above which the ordered 
phase disappears altogether. Two recent papers, by Nat- 
termann, Scheidl, Korshunov and Li (NSKL) |lq| , and 
by Cha and Fertig |ll| , cast doubt on the rccntrance pic- 
ture. The phase diagram they suggested is shown in Fig. 
1(b), w here the reentrance line TL(cr) disappears. NSKL 
p8| , po|j2l| further suggested that some sort of freezing 
phenomenon takes place below a certain temperature 

T*(ct) = 2aJ (1.5) 

[see the dashed line in Fig. 1(b)], which preempts the 
reentrance transition at T_ (er) < T* (a) found previously. 

The aim of the present paper is to expand the pio- 
neering ideas presented in Refs. Jl8| and Jl9| to unfold 
the physics which underlies the vortex-antivortex unbind- 
ing transition in the presence of the quenched disorder. 
There are two main extensions contained in this work as 
detailed below. 

First, we analyze quantitatively the equilibrium behav- 
ior of a single vortex in a background of quenched random 
dipoles. Analogy is made to two well-studied problems 
involving disorder: the random energy model p^ j, and 
a directed polymer on the Cayley tree p3fl . It is shown 
that the single- vortex problem has a glass transition at a 
temperature 

T g = J(W2) 1/2 , (1.6) 

below which entropy goes to zero, i.e., the vortex becomes 
localized at the lowest energy site. The free energy of the 
vortex is proportional to the logarithm of system size at 
all temperatures, with a prefactor which vanishes on the 
phase boundary shown in Fig. 1(b). 

Second, the dielectric and freezing properties of a di- 
lute gas of vortex-antivortex pairs (or molecules) are ex- 
amined in further detail, with particular emphasis on the 
spatial structure of equilibrium pair configurations. The 
freezing line T = T* in Fig. 1(b) is shown to be related to 
the loss of entropy of a pair over an area where the pair 



can be considered as isolated from other pairs of compa- 
rable size. If we fix the center position of the pair, the 
two vortices make up the pair freeze at T g . In the ordered 
phase, < T g due to the fact that the pair is allowed to 
explore an area much larger than its size and hence has 
a lower freezing temperature. Interestingly, freezing of 
pairs is not associated with a singularity in the free en- 
ergy of the system as a whole, and there is no real phase 
transition at T*. Disorder also generates random, zero- 
field polarization of the gas of pairs, which enhances the 
effective disorder seen by vortices separated by a large 
distance. This effect, which has been previously over- 
looked, shifts the horizontal phase boundary a = tt/8 in 
Fig. 1 (b) to smaller values of a . 

The outcome of these considerations can be turned into 
a set of renormalization group (RG) recursion relations 
which capture the average, large- distance properties of 
the system. Apart from some minor differences, the RG 
flow equations derived in this paper are in agreement 
with those of Ref. Q. To the extent that such a sim- 
plifying description offers a good approximation, a phase 
diagram of the form Fig. 1(b) is produced. The RG 
description is however not sensitive enough to rare fluc- 
tuations. The influence of rare fluctuations on some of 
the quantitative aspects of our results, such as the slope 
of the phase boundary as T tends to zero, remain to 
be studied. Qualitatively, though, the basic conclusions 
of the RG calculation are expected to be valid, as the 
modification of the bare interactions due to excitation of 
large-size pairs is relatively small in the entire ordered 
phase shown in Fig. 1(b). 

An interesting question is whether the system in the 
low temperature region above the <7 c -line has long-range 
glassy order. Our calculation of the dielectric suscepti- 
bility of a gas of pairs indicates that screening is present 
at all temperatures, despite localization in the orienta- 
tion of individual pairs below T g . This supports the idea 
that, in the disordered phase, vortex-vortex interaction 
at large distances is always short-ranged. Consequently, 
long-range glassy order in the phase field is not expected 
at any nonzero temperature due to finite energy cost to 
excite an additional vortex in the system. 

The paper is organized as follows. In Sec. II the 
Coulomb gas representation of vortices of the XY model 
is briefly reviewed. A qualitative discussion of vortex- 
antivortex unbinding is presented to highlight the out- 
standing issues. The problem of a single vortex interact- 
ing with quenched random dipoles is analyzed in Sec. III. 
Connection is made to the random energy model and to a 
directed polymer on the Cayley tree. In Sec. IV we exam- 
ine the behavior of a dilute gas of vortex-antivortex pairs 
of comparable size, under the influence of disorder. The 
calculation of the dielectric susceptibility and the zero- 
field polarization of such a gas is presented, as well as an 
analysis of fluctuations of pair density. A physical inter- 
pretation of the T» line is proposed. Section V contains a 
derivation of the RG recursion relations and results that 
follow from these equations. A discussion of the phase 
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diagram, singularity of the free energy, divergence of the 
correlation length, and the two-point phase-phase cor- 
relation function is presented. The main results of the 
paper are summarized in Sec. VI. Some of the technical 
aspects of the study are relegated to the four appendices 
at the end. 



(See Appendix A for more details on the derivation.) 
Here and elsewhere = r j — Tj is the displacement vec- 
tor between sites i and j, and ry = |r,j | is the distance. 
In addition to the usual core energy E c , a vortex interacts 
with a quenched random dipolc field q; = (a/27r)A(rj)xz 
through the potential 



II. COULOMB GAS FORMULATION 



A. Vortex Hamiltonian 



To set the stage, let us review brief ly th e steps leading 
to the Coulomb gas representation of ( |l . l[ ) . T he s tandard 
procedure is to take a continuum limit of (LI), which 
yields a quadratic Hamiltonian |^4|,|| , 

H= J - J d 2 r[V4>- a' 1 A(r)] 2 , (2.1) 

where a is the lattice constant. The two components of 
A are given by the disorder Aij on adjacent horizontal 
and vertical bonds, respectively. 

In the presence of vortices, the field cj>(r) is multi- 
valued. The vortex configuration is specified by a set 
of vortex charges rrii such that the phase advance along 
a closed path surrounding site i (or rather cell i) is given 
by 



= 2irmj 



In a system with periodic boundary conditions, neutral- 
ity J2 i mi = is satisfied. The gradient of the 4> field 
can now be decomposed into a rotation-free part and a 
divergence- free part, 



= X7<t> sw + ^m^z x (r - rj)/|r - r 
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(2.2) 



where (j) sw represents "spin- wave" fluctuations. The same 
procedure can be repeated for A, 



A = aW<f> + A r , 
where the potential 4>q satisfies 



aV 2 



= V- A. 



(2.3) 



(2.4) 



Inserting Eqs. (2^2) and (2J3) into ( |2.l|) , we obtain 
(apart from a constant) H — H sw + H v , where the spin- 
wave part is given by 



H sw = ^ J d 2 r(S7(j) sw ~ V0 O ) 2 
and the vortex part given by 



(2.5) 



H v = ^^(m 2 E c + rriiVi) — irJ ^""^ m^mj In — — . (2.6) 

i i^j 



V i = V(r i ) = 2TrjJ2<lr^/ri j . (2.7) 

From the above definition we have, in component form, 
(qi,a) = 0, {q%, a <li,p) = {a/2ir) 2 <j5rj8 a p. (2.8) 
Note that Vi vanishes when A is rotation-free. 

B. Pair- unbinding transition 

At sufficiently high core energies, at least, the gas of 
vortices in a charge-neutral system is expected to form 
one of the two phases described below. The first is a di- 
electric insulator, where ±1 charges bind to form pairs 
of charge-neutral molecules. This structure is low in the 
Coulomb energy, but also low in entropy due to binding. 
The second is a plasma with a finite density of unpaired 
(or free) vortices. This structure is high in the Coulomb 
energy but also high in entropy. In the absence of dis- 
order, both the Coulomb energy and entropy scale log- 
arithmically with distance in two dimensions. A simple 
energy-entropy argument then predicts a finite temper- 
ature transition for the unbinding of vortex-antivortcx 
pairs. This is also the temperature where the free energy 
of a single vortex goes to zero. An improved treatment, 
which takes into account reduction of the Coulomb en- 
ergy due to screening by other vortex-antivortex pairs, 
yields an exact description of the critical properties at 
the transition. In the plasma phase, there is complete 
screening of the Coulomb potential, so that interaction 
between distant charges become short-ranged. 

In the presence of quenched random dipoles, vortices 
may explore fluctuations in the disorder potential to 
lower their Coulomb energy, and hence become more nu- 
merous. This speaks for the reduced stability of the in- 
sulating phase. On the other hand, in the process of 
gaining potential energy, vortices become more localized, 
and this way loose entropy. The first insight one needs 
is how much energy a vortex can gain from the disorder 
by positioning itself at the right place. It turns out that 
this problem can be solved almost exactly, and the result 
again has logarithmic scaling with distance. The ampli- 
tude of energy gain from disorder is proportional to a 1 !' 2 
at low temperatures. Thus, when entropy is not a factor, 
excitation of free vortices is not expected below a certain 
critical strength of the disorder. 

As in the pure case, a complete treatment requires 
analysis of screening of the Coulomb potential due to 
other pairs of vortices present in the system. At high 
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temperatures, a pair is able to explore a large number 
of different disorder environment, which minimizes the 
difference between quenched and annealed disorder. The 
situation becomes different at low temperatures where, 
as in the random energy model, the equilibrium behavior 
of a pair is dominated by the lowest energy configuration 
in the area accessible to the pair. A crucial issue is thus 
to obtain the correct statistics of the pair when spatial 
and angular localization becomes important. 

With the above general picture in mind, we are in a 
position to perform the necessary calculations. 



III. SINGLE VORTEX 

In this section, we examine the behavior of a single 
vortex, confined in a box of linear dimension R 3> a. In 
the presence of disorder, the energy of the vortex depends 
on its position i, 



(/) = -T(lnz) = -c(T, s) In N + 0(ln In N), (3.7) 



Ei = E c + nJln(R/a) + Vi, 



(3.1) 



where Vi is given by ( |2 . 7| ) with the sum restricted to sites 
in the box. The variance and spatial correlations of Vi 
are given by 



(Vj 2 ) = 2it<tJ 2 ln(R/a) + 0(1), 
((Vi - Vj) 2 ) = 4naJ 2 hi(nj/a) + 0(1). 



(3.2) 
(3.3) 
(3.4) 



A simplifying approximation to the single- vortex prob- 
lem is obtained by setting the correlation of Vi to zero. 
The resulting problem is known as the random energy 
model (REM) |22|. It turns out that, for quantities of 
interest to us, correlations in the disorder potential only 
introduce minor corrections to the REM results. In the 
following we shall first discuss the REM and then an im- 
proved representation. 



A. Random energy approximation 



In the REM one considers the partition function 



JY 



= ^exp(-Xj/T), 



(3.5) 



where 



c(T, s) 



Here 



T + s/(2T), for T>T g (s); 
(2s) 1 / 2 , for T<T g (s). 



T g (s) = (s/2) 1 ' 2 



(3.8) 



(3.9) 



is the freezing temperature of the model. For T < T g , 
the entropy is no longer extensive in In A. 

The above result can be applied to the single vor- 
tex problem by substituting A — ► (R/ a) 2 , A — > 



2iraJ 2 ln(i?/a), and s — + iraJ 2 . From (3/7), we obtain 
the average free energy of the vortex, 



(F) 



E c + ttJ(1 - - In -, for T > T g ; 
irj 1 a 



E c + nj(l 



"°' N m -, for T<T q . 

it ' a 



(3.10) 



The corresponding freezing temperature is given by Eq. 
(1.6) (solid line in Fig. 2). The coefficient of the loga- 



rithm changes sign across the dashed line shown in Fig. 
2, which is precisely the phase boundary in Fig. 1(b) 
when renormalized values for J and a are used. (See 
discussion in Sec. V.) 

Below T g , the entropy of the vortex is no longer exten- 
sive in ln(R/a). In fact, it can be shown that only one or a 
few lowest e ner gy sites contribute significantly to the par- 
tition sum ( |3.5| ) in this regime (see Appendix B). Within 
the region bounded by the dashed line in Fig. 2, the typ- 
ical free energy of a vortex is positive, but there are rare 
realizations of disorder which give rise to a negative free 
energy. The probability for such events is a power-law 
function of R/a with a negative exponent. This fact is 
important when we consider pair excitations in Sec. IV. 



B. Correlations in the disorder potential 



where Xi, i = 1, . . . , N, are a set of random energy levels 
drawn independently from a gaussian distribution, 



ij}(x) = (27rA)- 1/2 exp(-a; 2 /2A). 



(3.6) 



The model has been analyzed in great detail by Derrida 
p2[ . Below we quote some of his results relevant for our 
discussion, and refer the reader to his original paper for 
further details. (See also Appendix B.) 

In the thermodynamic limit A —> oo while fixing the 
ratio s = A/ In AT, the average free energy is extensive in 
In A, 



The REM approach to the single- vortex problem is not 
completely satisfactory as it ignores spatial correlations 
in the energies Vi. This correlation has a simple origin 
(see Fig. 3). When we move the vortex from a site i to a 
site j, the change in the disorder potential is mainly due 
to a change in the local environment up to a distance of 
order nj, as contributions to Vi and Vj from quenched 
dipoles further away are nearly identical. This type of 
correlation can be easily coded using the Cayley tree, 
where each site is associated with a path on the tree. 
The potential on a site is made equal to the energy of 
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a path on the tree. Geometrical proximity is translated 
into hierarchical proximity on the tree. 

This representation can be made explicit using the fol- 
lowing construction, though details of it should be unim- 
portant for our conclusions. For any chosen site i, we 
divide the space into a set of rings of inner and outer 
radii R n -i and R„, respectively, such that a — Rq < 
R% < .. . < R rn ~ R, while keeping R n /R n -i = b con- 
stant. The potential at site i can be written as a one- 
dimensional sum, Vi = V i , where each term in the 
sum contains only contributions from dipoles within a 
given ring, i.e., 

"'-2-kj Yl q* •'<*/»■?*• 

Rn-i<r ik <Rn 



v 



(3.11) 



We now identify the nth ring with the nth node (branch- 
ing point) along the path i on the tree, where n increases 

from bottom to top. The energy of the node is given 

(n) 

by V t . Repeating the above procedure for a differ- 
ent site j, we obtain another sequence of energies Vj n ^ 
for nodes on the path j. The two paths join on level 
fiy = ]n(nj / a) / hxb. 

An intriguing fact about the random dipolar interac- 
tion is that the subsums constructed above are gaussian 
random variables with identical statistics, 



Thus all rings contribute equally to the sum Vi , indepen- 
dent of the radius of the ring. 

The Cayley tree problem discussed above has been an- 
alyzed in detail by Derrida and Spohn Q . Its properties 
are quite similar to the REM. In particular, the exten- 
sive part of the free energy is the same as in the REM, 
independent of the choice of b. In addition, moments of 
the partition functio n hav e the same dependence on N 
as indicated in Eq. (Bll), and the transition tempera- 
ture T n of the n-th moment is the same as in the REM. 
There are, however, differences in the amplitude of the 
ratio (z 11 ) / (z) n . This implies that the distribution of the 
f ree energ y, / — — TTnz, is not exactly given by Eqs. 
(B8) and (B9) for / significantly less than (/), but the 
difference should be small, as otherwise the behavior of 
(z n ) would be significantly different. 



IV. DILUTE GAS OF PAIRED VORTICES 

As mentioned in Sec. II. B, a quantitative study of 
the pair-unbinding transition must include a discussion 
of pair-excitations which modify the Coulomb interaction 
at large distances. This is usually done by employing 
a real-space RG procedure, to be explained in detail in 
Sec. V. A crucial step in the RG scheme is the calculation 
of the dielectric susceptibility and zero-field polarization 
of a gas of pairs in a certain size range, say between R 
and R + dR. This is the task to be carried out in this 
section. 



A. Lattice gas representation 



To treat a dilute gas of pairs of uniform size R, it is 
useful to separate the "internal" degrees of freedom of 
a pair, given by allowed configurations of the pair con- 
fined to a box of linear size R, from rigid translations of 
the pair over a distance greater than R. One way of im- 
plementing the idea is to impose a lattice with a lattice 
constant R. The lattice-gas representation is extremely 
handy owing to the following two properties of the sys- 
tem: (i) the disorder potential on a pair is essentially 
uncorrelated when the pair is translated over a distance 
larger than i?; (ii) interaction between pairs of similar 
size in the dilute limit can be approximated by a hard- 
core potential extending over a distance of the pair size 
R. These facts can be established following a similar line 
of reasoning as in the original paper by Kosterlitz and 
Thouless §. 

Let r + and r~ be the coordinates of +1 and —1 charges 
in a pair, respectively. The pair energy is given by 



E p = 2E C 



2TrJln{R/a) + V(r+)-V{r-), (4.1) 



(y/"))=0, (V} n) v} n,) ) =2n(jJ 2 (\Tib)5 n , nl . (3.12) where i?=| 



r | is the size of the pair. 

The rapid decay of correlations in the disorder poten- 
tial V(r + ) — V(r~) on a pair beyond a distance of order 
R comes from an observation made in Sec. III. B. The 
two charges which make up a pair interact separately 
with quenched random dipoles within a distance of or- 
der R from the pair center, but collectively as a dipole 
when more distant disorder is in question. Hence the 
random part of E p is dominated by disorder within a 
distance of order R from the pair center. (The remaining 
contribution from distant quenched random dipoles can 
be treated as a perturbation when necessary.) On the 
other hand, barring contributions from distant quenched 
dipoles, V(r + ) is quite independent from — V(r~) for two 
reasons. First, each potential is dominated by quenched 
dipoles in the immediate vicinity of the site in question 
(see discussion on the ring structure in Sec. III. B). Sec- 
ond, although the two charges are in the same disorder 
environment, when it comes to optimizing their (free) 
energies, they see opposite ends of the disorder energy 
distribution due to the difference in sign. Therefore, to a 
good approximation, we can repla ce E v by the sum of two 
single vortex energies of the form (3.1), each containing a 
random potential generated by quenched dipoles within 
a box of linear size R, independent from the other. 

The interaction between one pair and another is of the 
dipole-dipole form at large distances, which is small com- 
pared to E p and can be treated as a perturbation. The 
interaction becomes more complex when two pairs are at 
a distance R\ < R, but it is generally repulsive, with 
a strength of order 4nJ\n(R/Ri). (Note that the two 



5 



pairs should be arranged in such a way that it is not pos- 
sible to regroup them to form ±1 pairs of smaller sizes.) 
For simplicity, we shall replace the interaction by a hard- 
core potential of range R. In the dilute limit, the main 
effect of this interaction is to prevent more than one pair 
to take advantage of a particular favorable configuration 
(and the ones very close to it), which turns out to be a 
very important constraint at low temperatures [lq| . 

We are now in a position to define the lattice-gas rep- 
resentation. We divide the plane into a square lattice of 
cells, each of linear dimension R. Any given cell has at 
most one pair, and pairs in different cells do not interact 
with each other. The Boltzmann weight on an occupied 
cell can be written as y p z p , where 



Up = {R/a) 



-2irJ/T 



exp(-2£; c /r) 



(4.2) 



is the pair fugacity and z p is the configurational partition 
function of the pair attached to the cell. Since there is no 
interaction between different cells, the partition function 
of the system factorizes into a product of cell partition 
functions 1 + y p z p . In addition, average over all cells can 
be replaced by an average over the disorder, as each cell 
represents an independent realization. 

To apply the lattice-gas description to the system of 
pairs in a given size range, say between R and R + dR, 
we need to specify z p in more detail. For the discussion 
to be meaningful, dR should be small enough so that the 
pair fugacity y p can be regarded as a constant, but large 
enough so that individual charges in a pair are allowed 
to explore their own local disorder environment without 
been severely constrained by the specified range of pair 
size. Both criteria can be met by choosing dR ~ R. The 
configurational partition function of an occupied cell is 
given by 



E 



exp[- 



(r++: 
_R<|r+- 



")/2£cell 
■-\<R+dR. 



V(r+) - V(r~)- 
T - 



(4.3) 



The potential V(r + ) — V(r~) inside a cell has a spatial 
correlation of similar nature as the potential on a single 
vortex discussed in Sec. III. To simplify the calculation, 
we shall again make the random energy approximation 
where this correlation is ignored. The parameters of the 
REM applied to the problem of pairs are, 

N = 2n(R/a) 4 (dR/R), A = Ana J 2 \n(R/a). (4.4) 

For dR ~ i?, the freezing temperature T g for the pair in 
a cell is the same as the freezing temperature of a single 
vortex, Eq. ([L"q ). 



B. Pair density 

In equilibrium, the probability of finding a pair in a 
given cell is given by 



W = 



y p z p 



y p z p 



(4.5) 



For a dilute gas, the typical value of W is given by 
Wtyp ~ y p z p ,ty P , where z p ^ yp is the typical value of z p 
(see discussion in Ap pendix B). Combining Eqs. ( [3.7] ), 
~), jB6l ), and (|J), we obtain, 



a J R 



9i 



typ 



R\ -2*K{\-y/tofr) dR „ 

a J !T' f0rT<Ts 



(4.6) 



Here K = J/T. The exponent of the power-law changes 
sign on the dashed line in Fig. 2. 

Like z p , W has a broad distribution. Its mean value 
(W) deviates significantly from W typ for T < T g . Since 
the n-th moment of z p grows much faster than (z p ) n for 
sufficiently large n, it is not possibl e to calculate (W) 
by expanding the right-hand-side of (|4.5|) as a power se- 



ries of y p z p . Nevertheless, the average can be calculated 
by treating the cases y p z p < 1 and y p z p > I separately, 
as done in App endix C. R esult s of the calculation are 
given by Eqs. ( |C1 1 ) and ( C13 ) in respective tempera- 
ture regimes. For R 3> a, a power-law dependence of 
(W) on R is found, 



(W) 



Rs4-2irK+2n*K 2 dR 



R 



for T > T* 



#x4-7r/(2a) dR 



R 



(4.7) 



for T < T* 



The exponent freezes to a temperature- independent value 
below T*. 



C. Zero-field polarization 

The disorder environment in a given cell specifies a 
favorable configuration for a pair in the cell. The break- 
ing of rotational invariance thus yields a zero-field dipole 
moment, 



Po 



_ J2pexp[-E p /T] 



1 



y p z p 



(4.8) 



where p = r + 
sum in Eq. (4 



r~ is the dipole moment of the pair. The 
8) is restricted to t he internal degrees of 
freedom of the pair, as in Eq. (4.3). 

Due to statistical rotational symmetry, (po) = 0. Its 
variance can be calculated approximately from the fol- 
lowing consideration. Note that po is small when many 
distinct configurations contribute to the cell partition 
sum z p . It becomes large when the lowest energy config- 
uration (and nearby configurations with approximately 
the same orientation of p) dominates. Based on the dis- 
cussion of Appendix B, it is reasonable to assume that 



G 



the latter occurs whenever z p is significantly larger than 
i ts t ypical value, %,,typ- Replacing p inside the sum in 
(18) by the dipole moment of the ground state, we make 
an error with a probability of the order of Wtypj which is 
smaller than (W). This yields the estimate, 



with 



(|p | 2 )/i? 2 = (^ 2 )+O((W0 2 ). 



(4.9) 



The calculation presented at the end of Appendix C 
yields, for T < T», 



(W 2 ) ~ (1-T/T*)(W). 



(4.10) 



For T > T», (W 2 ) decays faster with R than (W). The 
distribution of |po| is expected to be broad. In particu- 
lar, for T <T g , where typically one or two configurations 
dominate the partition sum, the distribution of |po| is 
similar to the distribution of W . 

Let us now consider the correlation between po and 
the total dipole moment of disorder in the cell, 



J2 *• 



(4.11) 



iGccll 



Since po is mostly determined by the arrangement of the 
disorder in the immediate vicinity of the two charges 
making up the pair, we expect the contribution to po 
from q to be small, but the effect is important for later 
discussions. To estimate the contribution, let us consider 
a quantity po, which is the equivalent of po under the 
replacement q^ — » q^ = q^ — (a/R) 2 q. From the third ex- 
ample of Appendix A, we see that switching on q is equiv- 
alent to switching on a polarizing field E g = — 27r 2 Jq/ R 2 . 
Linear response theory then suggests, on average, a rela- 
tion of the form, 



Po - Po - 2tt xJq, 



(4.12) 



where \ is the average dielectric susceptibility of the gas 
of pairs, to be discussed below. 



X=(2T)-\(W)-(\p \ 2 )/R 2 ). (4.16) 
) we may rewrite the above 



Using Eqs. (4.5) and 
equation as, 



X 



1 r d(W) 



2T 



Vi>- 



dy P 



0({Wf 



Here we have used the identity 

y p d(W)/dy p = (W) - (W 2 ) 



(4.17) 



(4.18) 



The derivative in the above equation can be evaluated 
using Eq. ( |CTo| ) for T > T», and ( gig ) for T < T*. To 
leading order, the result reads, 



{2T)- 1 {W), for T>T* 
(2T,)- 1 (W), for T<T* 



(4.19) 



[Note that, in both cases, the coefficient in front of (W) is 
fixed by th e (eff ective) power-law dependence of (W) on 
y p . Hence (4.19) is more exact than what one might have 
expec ted from the approximate nature of Eqs. (CIO) and 

The dielectric susceptibility is finite down to T = 0. 
At T = 0, individual pairs can not respond to a weak 
applied field due to loss of entropy. The polarizability of 
the medium is a consequence of a finite density of states 
at zero pair energy. Pair configurations with a slightly 
positive energy in the absence of the field may acquire a 
negative energy if it is favored by the field, and hence be- 
come occupied. The opposite happens for the unfavored 
pair configurations opposing the field direction. 



E. The pair freezing temperature T» 



D. Induced polarization 

In the presence of a weak, constant external electric 
field E, a cell acquires an induced dipole moment due to 
pair excitation, 

_ Epex P [-(g p -p-E)/T] 
1 + L exp[-{E p - p ■ E)/T] 

The induced polarization Pj„d of the gas of pairs is given 
by the spatial average of Pind, or equivalently, the disor- 
der average, 



Pind = R 2 (pind) 
To the first order in E, we find 

Pind = XE, 



(4.14) 
(4.15) 



The change in the leading order behavior of the pair 
density p ~ (W)/R 2 at T* has a simple interpretation. 
Given the strong repulsive interaction between two pairs 
at a distance smaller than their size R, and the absence 
of correlation in the disorder potential on a pair beyond 
a distance of order R, it is reasonable to assume that 
clustering of pairs is rare in the dilute limit. The typi- 
cal distance between neighboring pairs is thus given by 
L = p- 1 ' 2 > R. Within an area of linear size L, we have 
typically one pair only. 

Let us first consider the equilibrium statistics of a sin- 
gle pair in a box of linear size L, taken to be arbitrary for 
the moment. The total number of configurations avail- 
able to the pair is N = (R / a) 2 (L / a) 2 , and the variance 
of the random potential, A = Ana J 2 \n(R/a). In the ran- 
dom energy approxim ation , the mean free energy of the 
pair follows from Eq. (|3.7D, 
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F p (L,T) 



2E C + 2itJ\a.{R/a) 

-2c(T, s)[\n{R/a) + ln(L/a)], (4.20) 



where 



2tto-J 2 \n(R/a) 
ln(R/a) + \n(L/a)' 



(4.21) 



For a fixed L, F p (L,T) increases with decreasing T, and 
locks to a constant for T < T g (s). At a fixed tempera- 
ture, F P (L,T) decreases with increasing L. 

The typical inter-pair distance L(T) is determined by 
the condition 



F p (L,T) = 0. 



(4.22) 



From the properties of F p mentioned above, we see that 
L(T) increases as T decreases, and locks to a constant 



L* for T < T*. Here T» 
consistently, with s* given by (4.21) 
result for T* agrees with (O). For T > T, 



is obtained self- 
at L = Z*. The 
we may use 



the high-temperature expression for c(T, s) in (4.20), and 
the condition ( 1.22| ) yields the following estimate for the 
number of pairs in an area of size R, 



{R/Lf ~ (R/a) 



4-2irK+2Tro-K 2 



exp(-2£ c /T), (4.23) 



in agreement with (Cll). The length L* satisfies 



(R/L,) 2 ~ (i?/a) 4 exp(-A/2T 2 ), 



(4.24) 



in rough agreement with ( |C12 ) for the number of pairs 
in an area of size R below T* . 

The physical meaning of the temperature T* is now 
clear. For T > T*, the entropy of a pair in a region of the 
size of inter-pair distance is finite and varies smoothly 
with T. This entropy is lost at T*. Therefore T* is as- 
sociated with the pair freezing. The length scale L*(R) 
is the smallest size of an area where one typically finds a 
negative ground state energy for a pair of size R. 

In contrast, the single- vortex glass temperature T g is 
associated with the lost of entropy for a pair when it is 
restricted to an area of pair size. [Note that (4.21) re- 
duces to the expression for a single vortex when we set 
L = R.] This temperature does not play a special role 
in the equilibrium behavior of a pair, where the relevant 
length scale is set by the inter-pair distance. Likewise, 
so far as the equilibrium properties of a dilute gas of 
pairs are concerned, the cell representation we employed 
is merely a convenient device for performing calculations. 

The equivalence of our results to those of Refs. |is|] 
and pl| implies that there is a simple connection be- 
tween the two approaches. In the work of NSKL and the 
more recent paper by Scheidl, calculation of thermal aver- 
ages were made under the "factorization-ansatz" , which 
assumes that pairs do not interact unless they take iden- 
tical positions. From the discussion of Sec. IV. A we see 
that the pair-pair repulsion extends to a distance of the 
order of pair size R. If there is no strong reason provided 



by disorder for clustering of pairs, the two approaches 
should differ only by a relative amount proportional to 
the pair density, i.e., the difference should show up only 
at order (W) 2 in the expressions for x, etc. This is pre- 
cisely what happens under the random-energy approxi- 
mation. In reality, due to correlations in the disorder po- 
tential, close to a very favorable configuration for a pair, 
there are other configurations which are nearly as favor- 
able, though pair-pair repulsion would forbid simultane- 
ous occupation of these configurations. The true density 
of pairs is thus expected to be somewhat smaller than 
the one calculated under the factorization ansatz or the 
random energy approximation. Nevertheless, from what 
we understand about the correlations, the qualitative be- 
havior of the system should be the same as predicted by 
the approximate calculations. In partic ular, no change in 
the exponent of the power-laws in Eq. (4/7) is expected. 



V. RECURSION RELATIONS AND RESULTS 
A. The RG transformation 

The knowledge we gained about a dilute gas of vortex- 
antivortex pairs can now be incorporated into a RG pro- 
cedure aimed at capturing the large-distance behavior of 
the Coulomb gas with disorder. This can be done explic- 
itly following an integration scheme used previously by 
Kosterlitz for treating the pure problem (|] . 

Consider a configuration {rrii} made up of two groups 
of charges. The first group, {mf}, consists of pairs of ±1 
vortices, each of size less than a cut-off size R. The sec- 
ond group, {mf}, consists of charges which do not fall 
into that category. (Note that our usage of the super- 
scripts "<" and ">" is the opposite of the one familiar 
in m ome ntum-space RG.) The total energy of the system, 
Eq. (2.6), can be rewritten as, 



#v(W) 



H v ({m<}) + H v ({m>}) 
+H int ({mf},{mf}), 



where 



hit 



E>(r n ) 



(5.1) 



(5.2) 



describes the interaction between the two groups. Here 
p„ is the dipole moment of pair n in the first group, r n 
is the center position of the pair, and 



E>(r) = 2njJ2> 



r - r< 



(5.3) 



is the electric field at r due to the presently unpaired 
charges in the second group. 

The partition sum over the paired charges is given by 

E< = cxp{-[H v ({m<}) + H int ]/T). (5.4) 
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Writing E*- = ^ exp(—SH/T), where Sq is the parti- 
tion function at iJj nt = 0, we obtain, 

SH({m>}) = -T]n{exp(-H int /T)) , (5.5) 

where (-)o denotes thermal average with respect to 
H v ({mf}) . Treating H lnt as a perturbation, we can write 
5H in a more suggestive form, 



SH = - d 2 r[P (r) 



ip ind (r)].E>(r)+0(|E>| 3 ) 



Here Po(r) = P"^( r — r n))o is the zero-field polar- 
ization of the paired charges in the absence of the in- 
teraction term H lnil and Pi, lc i = P Po is the induced 
polarization of the paired charges due to the field E>. 
[Note that P(r) is defined in the same way as Po(r) ex- 
cept that the thermal averaging is taken with respect to 
H v ({m<}) + H int .} 

The renormalization group idea is to take the cut-off 
size R as a running parameter, and perform the elimina- 
tion of paired charges {mf} in a step by step manner, 
so that each time one needs to deal with pairs in a nar- 
row size range R to R + dR only. The necessary calcu- 
latio ns ha ve alr eady been done in Sec. IV. Substituting 
Eq. ( 4.15 ) into (5J3), we obtain, 



SH 



d 2 r[Po-E> + -x|E>| 2 ]. 



(5.7) 



This is nothing but the field integral version of the 
Coulomb energy (|2.6|), and hence can be incorporated 



into H v ({mf}) by redefining the parameters J and er of 
the model. 

The change in J can be obtained with the help of the 
first example in Appendix A. One thing to note is that 
the integral over |E > | 2 in (5.7) exclude regions of size R 



around each mf charge, since paired vortices should not 
be found in these areas. This leads to the identification 
b = R in Eq. (AS). The new effective parameters are 
given by, 



E C ^E C = E C + 47r 3 xJJln(i?/a), 



and 



J -> J- 



J- 



■ 4tt 2 x, 



(5.8) 



(5.9) 



equivalent to Eq. ( A4) |25|]. The extra term in E c merely 
accounts for the fact that screening from this group of 
pairs is effective only at distances larger than R. 

In Sec. IV. C, contribution to the zero-field dipole mo- 
ment po of a cell due to disorder within the cell was 
calculated. More distant disorder contributes to po by 
acting as an additional polarizing field . W hen the latter 
contribution is substituted into Eq. (5.7), we see that, 
with the help of the second example in Appendix A, the 
interaction strength between {mf } and quenched dipolcs 



is reduced by a factor 1 — 2ir 2 x J- Combining the zero- 
field polarization P = Po/R 2 of pairs with the disorder 
polarization Q = q/i? 2 , we obtain the effective disorder 
that couples linearly to E> in the Hamiltonian H v + SH 
for the {mf } charges, 



Q cff = (1-47T 2 V/)Q + P , 



(5.10) 



where Po = 
around Eq. 




is independent of Q [see discussion 



Equation ( 5.10] ) shows that pair excitations modify the 
(5 6) quenched disorder seen by {mf} charges in two different 
ways. The first effect is the screening of the interaction at 
distances larger than the pair size R, w hich can be taken 
into account by a redefinition of J, Eq. (5.9). The second 
effect is the generation of additional disorder. Since Po 
is independent of Q, we obtain an additive contribution 
to the variance of disorder, a. Writing JQ = </Q c ff, and 
using (|Q| 2 ) = a/(2n 2 R 2 ), we get 



+ 2ir 2 (W 2 ) +0{(W} 2 ). 



(5.11) 



In deriving the above expression we used the fact that 
the difference between the variance of Po and that of Po 
is of order (IV) 2 . Using the result for (IV 2 ), we see that 
the change in a is proportional to (W) for T < T„, but 
of higher orde r fo r T > T* . 



Equations (5.E) and (5.11) can be expressed in the 

For con- 



ae 



usual differential form by writing R 
venience, we introduce a dimensionless quantity Y(R), 
such that 2irY 2 dR/ R = (IV) gives the number of vortex- 
antivortex pairs of size between R and R + dR, in an 
area of size R 2 and averaged over the whole system. For 
T > T* or K- 1 = T/J > K~ l = 2a, we have, 



dK^/dl = 4tt 3 V 2 , 
dcr/dl = 0, 
dY/dl = (2-tiK 



iraK 2 )Y. 



For T < T* or R- 1 < R-\ we have, 

dK~ l /dl = 2n 3 a- 1 K- 1 Y 2 , 

da/dl = 2tt 3 (2 - a^R-^Y 2 , 

dY/dl=(2-^-)Y. 



(5.12a) 
(5.12b) 
(5.12c) 



(5.13a) 
(5.13b) 

(5.13c) 



Note that, in writing the above equations, we only kept 
terms up to order V 2 . The flow equations for Y follow 
from the power-law depende nce o f the pair density on 
pair size as given by Eqs. ( |Cll| ) and (C13 ). A term 
of order l/l inside the brackets in ( 5.13c ) has been ne- 
glected. 

A few remarks concerning the a bove recursion rela- 
tions are in order. For T > Eqs. (5.12) i s iden tical to 
those of previous authors |||l2| . From Eq. ( 5.11 ) we see 
that there is a renormalization of a even in this regime, 
but the effect is of higher order than V 2 . The change 
of the flow equations for T < T* was pointed out ear- 
lier by NSKL |lq|. The renormalization of a, though not 
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recorded previously, has also been obtained by Scheidl 
using a different approach Q . 

In the absence of disorder, Y is equal to the "rescaled" 
single vortex fugacity, {R/a) 2 ~' Kj / T exp{— E c /T), in the 
dilute limit. When disorder is present, relation between 
Y and the core energy E c is more compl icate d. The bare 
value of Y can be obtained from Eqs. (|Clll) and jC13| ) 
for T > T* and T < T* , respectively. It has a finite limit 
Yb at T = 0. For small a, Y 2 ~ er 1 / 2 exp(— c/cr), where c 
is a positive, model-dependent number. At small values 
of T, the bare value of Y increases from Y by an amount 
proportional to T 2 . 



B. Phase diagram and thermodynamic properties 

1. Constants of RG flow 

Seemingly complex at first sight, the flow equations 
( |5.12 ) and ( 5.13Q have in fact the same structure as their 
a = counterpart [|| . The fixed points of the flow are 
located on the Y — plane in the three-dimensional (3D) 
parameter space spanned by A" -1 , er, and Y. They are 
stable in the region enclosed by the dashed line in Fig. 
4, but unstable outside the region. Points on the dashed 
line are hyperbolic fixed points which describe the pair- 
unbinding transition. 

It turns out that the flow equations are completely in- 
tegrable. In the region T > [to the right of the dotted 
line in Fig. 4], 



a = const. 



(5.14a) 



is obviously a constant of the flow. On the T < T* side, 
the corresponding first integral is given by 



A — aK = const. 



(5.14b) 



as can be easily verified using Eqs. (5.13a) and (5.13b). 
These "streamlines" of the flow are illustrated in Fig. 4. 
The second constant of the flow is given by 

Y 2 -n^iR- 1 - |erA + | In A) = const. (5.15a) 



for T > T*, and 



Y 2 -^{K- 1 



a H — In A) = const. 
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(5.15b) 



for T < T». A unique trajectory in the 3D parameter 
space is specified when one combines (5.15) with (5.14), 
with the constants fixed by bare values of the parameters 
involved. 



2. Phase 



The original XY model has only two parameters, 
Kg 1 = T/J and erg. The bare value of Y, Yb, is a func- 
tion of Kb and <jb- The phase boundary of the model 



is determined by the condition that the RG flow ends on 
the dashed line in Fig. 4. Since the flow takes both a and 
A -1 to larger values, this phase boundary lies within the 
area bounded by the dashed line, as illustrated in Fig. 5. 

The question of reentrance of the disordered phase is 
whether the upper-left part of the phase boundary illus- 
trated in Fig. 5 contains a piece with a positive slope. 
Although earlier calculations which led to the prediction 
of a reentrance transition are not to be trusted, it is ac- 
tually difficult to rule out s uch a possibi lity fr om the 
new RG flow equations ( 5.13| ). From Eq. ( 5.13b| ) we see 
that, at a given Y 2 , the increase in the effective disor- 
der becomes slower as temperature increases. Hence for 
a fixed Yb, the bare value of a which flows to the fixed 
point value a — 7r/8 increases with increasing T/J. On 
the other hand, in the original XY model, the bare pair 
density Yj at a fixed a is expected to increase with tem- 
perature, too. The calculation presented in this paper 
is not quantitative enough to assess the two competing 
effects to reach a precise conclusion on the shape of the 
low temperature part of the phase boundary. 



3. Approach to criticality 

The critical behavior around the transition is con- 
trolled by the RG flow close to the relevant hyperbolic 
fixed point. As an example, let us consider flows along 
a particul ar co ntour in Fig. 4. Substituting Eqs. (5.14) 
into Eqs. fl5.15| ), we get a flow pattern depicted in Fig. 6. 
The curve c onsists of two pieces, one from Eqs. ( 5.14b| ) 
and (5.15b ) for A -1 < A^ 1 , and the other from Eqs . 
Q 5.14aD and (pU5a|) for A" 1 



> K~ x . Since Eqs. (jfUj) 
do not involve Y z , any vertical translation of the curve 
shown in Fig. 6 is also an invariant of the RG flow. This 
family of curves can be parametrized by the minimum 
value of Y 2 on each curve, D. 

In the ordered phase and at the transition, the scaled 
pair density Y 2 eventually reaches zero on large length 
scales. The corresponding RG flow follows a curve with 
D < 0, with the marginal case D = reserved for the 
transition. On the disordered side but close to the tran- 
sition, the RG flow follows a curve with a small positive 
D. Along such a curve, the value of Y 2 first decreases 
as the running parameter I increases, becomes almost 
stationary around the minimum of the curve, and then 
grows rapidly to large values beyond I = l + . The length 
£ + = aexp(Z + ) sets the scale where the correlation func- 
tion (see Sec. V. C below) turns from a power-law to an 
exponential decay. 

The dependence of on D can be obtained by exam- 
ining the flow in the vicinity of the minimum of the curve 
at A" 1 = KJ 1 (see Fig. 6). Let X = A" 1 - KJ 1 be 
the distance from the minimum. Generically, for small 
A, we have 



Y = XX' 



D, 



(5.16) 
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where A is proportional to the curvature at the minimum. 
It is easy to check that, in this case, t he coe fficie nt of Y 
on the right-hand-side of either Eq. (5.12c) or (5.13c), 
which ever applies, is linearly proportional to X. Using 
(5.16), we can write the flow equation for Y as, 



dY/dl^j-^Y 2 -D) 1/2 Y, 



(5.17) 



where 7 is a nonuniversal number which depends on the 
location of the fixed point. The parameters D and 7 can 
be scaled away using the substitution Y — > D X / 2 Y and 
I — > 7_D -1 / 2 Z. This yields a correlation length 



= aexp(^ + ) ~ aexp(7/L> 1/2 ). 



(5.18) 



(I 



The dependence of D on the shift of the bare parame- 
t crs f rom t heir c ritical values can be obtained by solving 
~14|) and fl5.15p . 

There is, however, a special case where the curvature 
of the curve at the minimum vanishes, invalidating the 
above analysis. This happens at the point S (and cor- 
respondingly S') in Fig. 5, where the T*-line meets the 
phase boundary. The minimum of the curve shown in 
Fig. 6 is now at the meeting point of the high and low 
temperature segments. This is also the inflection point 
of each of the two curves. The function around the min- 
imum takes the form, 



Y 2 



WX\ 



D, 



(5.19) 



dF v 
d\nR 



-2ttTY 2 /R 2 , for T>T*; 
-2irT*Y 2 /R 2 , for T < T». 



(5.24) 



The total vortex contribution F v to the free energy 
density of the XY model is obtained by integrating dF v 
over R. Let F VyC be t he vor tex free energy density at the 
transition. Equation ( 5.24 ) then yields, 



F v -F v 



R~ 3 dR(Y 2 - Y 2 ), 



(5.25) 



where Y C (R) is the value of Y(R) at the transition. A full 
analysis of the integral is quite involved, but the follow- 
ing consideration should yield a correct estimate of the 
expected singular behavior. 

Approaching the transition from the ordered phase, Y 
flows to at R = 00. The difference Y 2 (R) - Y 2 (R) is 
expected to remain constant, say equal to D < 0, up to 
R = and then goes to 0. Truncating the integral at 
i? = £_, we obtain, 



F v ._ - F v 



-Da' 2 



DC 2 - 



(5.26) 



The crossover length £_ has the same behavior as £,+ 
except that one should replace D by — D in Eqs. (5. IS 



and (5.21), as the case may be. From the disordered side, 
the story is the same up to R = , but beyond that Y 2 
becomes of order 1. Hence we expect 



where A is another constant. The flow rate of Y is now 
proportional to X 2 Y . Consequently, we have, 



F v 



F v 



-Da- 2 - C 



(5.27) 



dY/dl = 7 _1 (y 2 — D) 2 / 3 Y. 



(5.20) 



where 7 is yet another number. The proper scaling in this 
case is Y — ► D l / 2 Y and I — > jD~ 2 / 3 l. A new dependence 
of the correlation length on D follows, 

~ aexp( 7 /£> 2 / 3 ). (5.21) 



4- Free energy 

Let us now turn to the behavior of the free energy in 
the vicinity of the transition. Following the discussion 
of Sec. IV. A, we may write the contribution to the free 
energy per unit area from pairs in the size range R to 
R + dR as, 



dF v ~ -T(\n{l + y p z p ))/R 2 
Comparing the above equation with 
d{dF v ) 



Vp- 



= -T(W}/R 2 . 



(5.22) 
we obtain, 

(5.23) 



By making analogy to Eqs. (4.17) and (4.19), and using 
the definition of Y, we find, 



The singularity of F y ± are thus related to the singular- 
behavior of £±. In both cases, it is an essential singular- 
ity. 

Finally, it is interesting to see if the pair-freezing line 
T* in the ordered phase corresponds to another singular- 
ity of the vortex free energy. The analysis presented in 
Appendix D suggests that this is not the case. We should 
however note that the absence of a true glass transition 
in our model is a result of the very special type of func- 
tional dependence of relevant quantities on temperature 
and pair size, and hence might be susceptible to various 
types of perturbations. 



C. Two-point correlation function 

Consider now the two-point phase-phase correlation 
function, 



C(rij) = (exp[z(<£i - fa)]), 



(5.28) 



where the average is taken over thermal and then disor- 
der fluctuations. When only spin- wave contributions to 
<f> are taken into account, one obtains a po wer- law decay 
of C with distance r,-j, as shown in Eq. (1.3). Vortex- 
pair excitations "soften" the spin-waves, and lead to a 
faster decay of the correlation function. On grounds of 
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renormalizability of the model, one expects that, in the 
ordered phase, 



C(r) 



(5.29) 



at large distances, but the exponent rj differs from r] sw 
in that the renormalized values and Coo should be 
used, 



n 



1 

27 



{Kj + (Too). 



(5.30) 



The above conjecture for the average behavior of C(r) 
can be verified by a direct calculation. Since the spin- 
wave and vortex fluctuations decouple, and the disorder 
which enters H sw is orthogonal to the disorder in H v (see 
Sec. II. A), we may write, 

C(r) = C sw (r)C v (r), (5.31) 

where C v (r) is the correlation of phases <p v generated by 
vortex excitations, 

Cv(r-ii) = (exp[z(0 ViJ - v ,i)])- ( 5 - 32 ) 

Followi ng K osterlitz ||], we calculate the right-hand- 
side of ( 5.32 ) by successive elimination of vortex- 
antivortex pairs, starting from the smallest pair size. The 
formula below gives the difference in the phase at two 
sites i and j, generated by a single vortex- antivortex pair 
at r , 



(27rJ)- 1 (pxz)-E(r p ), 



(5.33) 



where p is the dipole moment of the pair (assumed to 
have a magnitude R much smaller than r) and E(r p ) is 
the electric field at r p due to a +1 charge placed at site 
i and a —1 charge placed at site j. The variance of the 
phase difference generated by all pairs in the size range 
R to R + dR can now be e asily calc ulate d with the help 
of the equivalence of Eqs. (]A6|) and (A8), 



<(0v 



2tt(W) Hnj/R). 



(5.34) 



Going from R to R + dR, C v (r) is reduced by a factor 
cxp[-i(W> v , - VJ ) 2 M - (r/R)- dv . (5.35) 



Comp arin g the c oefficie nt of the logarithm in ( 5.34| ) with 
Eqs. ([l9|), (|t|), and ([Til]), we obtain, 



da), 



(5.36) 



which is nothing but the differential form of ( |l.4| ). This 
confirms the i ntuiti ve idea that the asymptotic value of 
r\ is given by ( |5.30 ). On the transition line, r\ decreases 
monotonically from the value 1/4 for the pure case to 
1/16 at the zero-temperature transition point p8[ . 

Inspecting Eq. (5.33), we see that spatial and angu- 
lar localization of the vortex-antivortex pairs may lead to 



rare, but significant deviation of the thermally averaged 
correlation function at two fixed sites i and j, 



C v (ri,rj) = (exp[i(0 v ,i - <t> v ,j)]) 



thermal 1 



(5.37) 



from its average value, C v (ry ). Such a fluctuation comes 
about when we are looking at large-size pairs which have 
strong density fluctuations (after thermal averaging) on 
the scale r^, and that the presence or absence of such a 
pair in the region surrounding i and j will make a signifi- 
cant change to the phase difference (p v ^ — <fi v ,j- Below T*, 
the typical density of these pairs is significantly smaller 
than the average density. Consequently, the typical value 
of the exponent, ??typ> can be somewhat smaller than its 
average value r/. 



VI. SUMMARY 

The main conclusions of the present work can be sum- 
marized as follows. When the variance a of the random 
phase shifts is smaller than a critical value <r c , the quasi- 
long-range order of the 2D XY model survives at suffi- 
ciently low temperatures. The value of a c is nonuniversal 
but should be smaller than n/8. The two-point phase- 
phase correlation decays algebraically with distance in 
the entire ordered phase up to the transition. At the 
transition, the exponent of the power-law decay lies in 
the range between 1/16 and 1/4. Approaching the tran- 
sition from the disordered side, the correlation length di- 
verges exponentially with a — i- or — |-power of the dis- 
tance from the transition point. The free energy exhibits 
an essential singularity on the phase boundary. The low 
temperature region of the disordered phase at a > er c is 
not expected to show glassy long-range order. 

The behavior of vortex-antivortex pairs inside the or- 
dered phase is quite interesting. In contrast to the pure 
case, there is a finite density of such pairs at zero tem- 
perature. As the temperature T increases, the pair den- 
sity increases initially slowly up to T — T„, and then 
grows rapidly as entropy comes into play. Following an 
opposite sequence, as T is lowered from the transition 
temperature, large-size pairs undergo various degrees of 
localization both in space and in its angular distribution. 
However, in the Coulomb gas language, a finite suscep- 
tibility for the gas of pairs is found at all temperatures. 
Localization also introduces a zero-field random polariza- 
tion of the gas of pairs, which has the effect of enhancing 
disorder seen by large-size pairs. 

Much of the qualitative aspects of our results agree 
with those of NSKL flf] and of Cha and Fertig g9), 
though there are minor but important differences on a 
quantitative level. Technically, the analogy introduced 
here to the random energy model offers a simple way 
to understand some of the subtle features introduced by 
the disorder potential, responsible for the failure of previ- 
ous calculations based on a small-fugacity expansion. A 
somewhat surprising result is that no singularity in the 
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free energy is found at T*. In fact, we have shown that 
an expansion with respect to the pair-fugacity y p fails 
for any a > 0. The whole analytic structure of the free 
energy as a function of T and a remains to be explored. 

A new quantity which appeared in this paper is the 
single- vortex glass temperature T g . This temperature 
also signals localization in the angular distribution of a 
pair when translation of the pair over a distance larger 
than its size is forbidden. As we have seen, this temper- 
ature plays no special role in the thermodynamics of the 
XY model. Nevertheless, one might contemplate possi- 
ble changes of dynamical behavior at T gi an issue to be 
studied further. 
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APPENDIX A: ELECTROSTATICS IN TWO 
DIMENSIONS 

In this Appendix we collect some useful formula from 
electrostatics in two dimensions. We adopt the conven- 
tion that the unit of charge is 1. The electric potential 
due to a +1 charge at origin is given by 



V(r) = -e _1 ln(|r|/a) 



(Al) 



where e = (2irJ) 1 is the bare dielectric constant. The 
electric field E = — W satisfies 



V • E = 27re _ V( r ) 



(A2) 



where p(r) is the charge density at r. 

In a dielectric medium, the induced polarization P = 
X E, where \ is known as the dielectric susceptibility. 
The displacement vector D = eE + 27rP satisfies 



V D = 27TP/0), 



(A3) 



where p/(r) is the density of free charges at r. Writing 
D = e E, we obtain 



2tt X . 



(A4) 



Let us now consider three examples encountered in the 
main text. In the first example we establish the equiva- 
lence of two expressions for the energy of a charge-neutral 
system. The electric field generated by a set of point 
charges mt located at r t is given by 



E(r) = 2vrJ^\ 



Consider now the integral 
1 



£ 



8ir 2 J 



d 2 r IE! 2 . 



(A5) 



(A6) 



Writing 



|E| 2 = W • E - V • (VE), 



(A7) 
2|) , we obtain, 



and using the Gauss's theorem and Eq. 



£ = ^2 ™> 2 E C - nJ^rriimj ln(|r j; - r 3 \/b). (A8) 

i ijtj 



Here 



2 



d 2 r\v\- 2 



r\<b 



(A9) 



is t he "core energy" of a unit charge. [The divergence 
of ( |AS| ) at small distances is cut off by the existence of 
a lattice.] Note that, in continuum, b is an arb itrary 
parameter which does not influence the result ( AS ) . 

As the second example we consider an expression for 
the interaction energy between a set of charges to, at 
and quenched dipoles q,, at Uj. Let E(r) be the electric 
field generated by all the charges, and E<j(r) be the field 
due to the dipoles but excluding those within a distance 
R from r. Consider now the integral 



Si 



1 



d 2 r E • E rf . 



(A10) 



4tt 2 J, 

We now write Ed as a sum over the field due to individ- 



ual dipoles. Using again (A7) and the Gauss's theorem 
separately for each j, we obtain, 



Si = ^ m i V d (r. 1 ) + ^Uj, 



(All) 



Here Vd(r) is the potential due to dipoles outside a circle 
of radius R centered at r, and 



1 2 



1 f 2 * 

— dej) R- E(R,- + R 
!7r Jo 



q 3 R 

R 2 ' 



(A12) 



where R = i?(cos <f>, sin 4>) . Elementary calculation yields 
Uj = |qj : • E(Rj). It follows that the second sum on 



the right-hand-side of ( All ) is — | times the first sum. 
Hence 



Si = -^miVd(rj). 



(A13) 



Our final example concerns the field Eq inside a circle 
of radius R generated by a medium with a permanent 
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polarization Q which fills the circle. To by-pass an ex- 
plicit calculation, we use the result that, in a uniform 
external field Eo, the field inside such a circle filled with 
a polarizable medium of dielectric constant ei is given by 



Ei 



2c 



ei 



-Eo, 



(A14) 



where e is the dielectric constant of the medium out- 
side the circle. The polarization inside the circle is 
Q = (ei — e)E 1 /(27r), and the the field it produces is 
Eq = Ei — Eq. Using the above results, we obtain, 



E Q = -(tt/c)Q = -2vr 2 JQ, 
which is uniform inside the circle. 



(A15) 



APPENDIX B: FREEZING AND LONG- TAILS IN 
THE RANDOM ENERGY MODEL 

As shown by Derrida [M , the freezing transition in the 
random energy model can be understood as a switchi ng o f 
terms which contribute most to the partition sum (3.5). 
For T < T g , the lowest of the N energies Xi dominates 
z, while for T > T g , typically a finite fraction of the N 
energies contribute significantly to z. (The word "typi- 
cal" refers to events which occur with a large probability, 
and "rare" refers to events with a very small probability.) 
This can be seen more explicitly as follows. 

In a given realization of the disorder, the random 



energies fall into a band 



^ Xi Xr\ 



where 



%mm = noin^-fxi} and a; max = maxj{a;j}. Introducing 
the integrated density of state, Af(x), which gives the 
number of levels with Xi < x, one can rewrite Eq. ( |3.5| ) 
as 



POO 

(-x min /T)+ / exp(~x/T)aW. (Bl) 



z = exp 



The contribution from the lowest energy level has been 
isolated from the rest. When the total number of lev- 
els is large, there are typically many levels in an interval 
dx ~ T, so that one may replace dJ\f by its mean, 



dAf ~ N^{x)dx. 



(B2) 



Equation (B2) fails when Af(x) < 1. This happens 
for x < xq, where xq is d eterm ined by the condition 
Af(xo) — 1. From Eqs. (B2) and (|3.6|) we get 



(B3) 



xq = -(2s) 1/2 lnAr + 0(lnlniV). 



Substituting (B2) into (Bl), and restricting the integral 
to x > xq, we obtain 



exp(-a; min /T) + z typ 



(B4) 



where 



^typ = N [ 



dxip(x) exp(— x/T). 



(B5) 



It is easy to show that the typical value of x m i n is given 
by xq , while the typical value of the partition function is 
given by 



z typ ~ exp(-(/)/T). 



(B6) 



Thus the contribution from the lowest energy level to z 
becomes significant in a typical realization of the disorder 
only when T <T g . 

On the other hand, fluctuations of z far away from Zt yp 
(say z > 2z typ ) are dominated by fluctuations of x m i n 
at all temperatures. This is especially so for T < T g , 
where the fluctuations of z are typically of order z typ . To 
characterize this behavior more precisely, let us consider 
the distribution of x m - ul , denoted by ipmin(x). For the 
minimum of the N energies to be greater than a certain 
number x, all N energies must be greater than x. Hence 
we have 



ip(y)dy] 



N 



(B7) 



For x <§; — v A, the integral on the right-hand-side is very 
close to one, in which case one can write, 



d f x 
ipmin(x) = -— cxp[-iV J ip{y)dy]. 



(B8) 



For x < xq , the argument of the exponential is less than 
1, and hence 



Ipmin(x) ~ Nlp(x) 



exp 



)lnJV . (B9) 



The distribution of x m - ln gives us an idea about the 
high-end of the distribution of the partition function z, 
where we can write z ~ exp(— x m - m /T). On a log-log plot, 
the local slope of the distribution is essentially given by 



CM 



din 
din; 



T T\nz 

Tn X 



(BIO) 



which is a slow- varying function of z. This implies that 
the distribution of z has a long tail. High moments of z 
are sensitive to the tail of the distribution. Using 
we find, 



(z)> 
where 



1, for T>T n ; 

Ar (n-i)[(r n /T) 2 -i ]) forT<T?i 



(Bll) 



(B12) 



The result agr ees with an exact calculation by Derrida 
starting from (|3~5|). 
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APPENDIX C: MEAN PAIR DENSITY AND 
FLUCTUATIONS 

The disorder average of the cell occ upa tion number W 
can be calculated by expanding Eq. (4.5) as a power se- 
ries of ypZp for y p Zp < 1, and a power series of (y p z p ) _1 
for ypZp > 1. Denoting by P(z) the probability distribu- 
tion of z p , we obtain, 

oo oo 

(w) = j^-i)- 1 /! w + £(-ir h(n), (ci) 



n=l 



n=0 



where 



/.I'D = //,',' / ^ z n P(z)dz, 
la 



h{n) = y~ 



Z - n P(z)dz. 



(C2) 
(C3) 



1 /Vp 



In principle, evaluation of the integrals I\ and I2 re- 
quires full knowledge of P(z), which we do not have at 
hand. On the other hand, as we discussed in Appendix B, 
the large- z tail of P{z) is due to fluctuations of the min- 
imum energy x m i n . Thus for large z the substitution 



P{z)dz — > ip m i n (x)dx ~ Nip(x)dx, 



(C4) 



with z = exp(— x/T), yields a good approximation. The 
integral I2 is now readily calculated. The result reads, 



I 2 (n)~—y p exp(^ 




(C5) 



where erf(u) = 27r x / 2 J Q " exp(— x 2 )dx is the error func- 
tion. Here 



T* = -A/(Tlny„) 



which coincides with (|l.5| ) in the limit i?, — > oc. 
The integral 7i (n) can be written as 

h(n)=y;(z n )-I 2 (-n). 



(C6) 



(C7) 



To obtain its leading order behavior, we need to examine 
which part of the distribution P{z) contributes most to 
the average (z n ). For T > T n = n l ' 2 T g , the main contri- 
bution comes from the central part of P(z) around z p ,typ, 
so that 



h(n) ^ (y P (z p )) 



(C8) 



For T < T n and T < nT*, the main contribution comes 
from the tail of P(z) at z > y" 1 . In this case, the two 
terms on the right-hand-side of ( |C7| ) almost cancel each 
other. The main contribution to I\ thus comes from P(z) 
around z = y^ 1 , where again the approximate expression 
for the tail of P(z) can be used. This yields 



h{n) ~ —y p exp(^ / 



1-erf 



(< 



/ 



'A' 
~2, 



(C9) 



When T n < nT*, the leading order behavior switches 
from (C8) to (|C9| ) at T„, though (09) appears as a sub- 
leading order term in the temperature range T n < T < 
nT*. (This is due to the fact that, around T = T„, 
(z n ) picks up significant contributions from z ~ z p.typ 
and from the tail at z > Up 1 -) On the other hand, 
when T n > nT,, there is an intermediate temperature 
range nT* < T < T n where 1\ is dominated by con- 
tributions from z between z P ,typ and Up 1 - In this case, 
I\(n) ~ y™(z n ), which esse ntial ly coincides with the ex- 



pression (C9). For n = 1, (C9) is valid at all tempera- 
tures. 

A careful analysis of the above expressions is quite 
cumbersome, but the following observations are useful 
and suffice for our purpose. 

(i) Leading order behavior. For T > T*, (W) is domi- 
nated by Ji(l), 



(W) ~iVy p exp(A/2T 2 ). 
In terms of parameters of the model, we have 
2E, 



(CIO) 



f}\ 4-2-rK+2tt<jK- 1 rlR 

4- ( cu ) 

a J R 



For T < T», I\(n) and hin) all contribute. Using the 
asymptotic expression erf(u) ~ 1 — 7r _1 / 2 u _1 exp(— u 2 ) 
at large u, we obtain 



< w > c " , (e)b^-^- a ^- 



(2ttA) 

or more explicitly, 

T\ / 2cr 



Here 



V2 / i?x4-7r/(2a) ^ 



\rj \ln(i?/a)/ la/ i? 



E 



(-1)™ 7TU 

n + u sin(7ro) 



(C12) 



(C13) 



(C14) 



The crossover from ( plO ) to ( C12j ) is not sharp, but oc- 
curs over a temperature range of order 



ST ~T4<r/]n{R/a)] 1/2 . 



(C15) 



[The apparent divergence of A(T/T*) at T = T* can be 
removed by separating out the contribution from Ji(l). 
This procedure yields a full description of the crossover, 
which we shall not elaborate here.] Note also that, since 
A(0) = 1, (W) has a finite value at T — 0, proportional 
to the density of cells with a negative ground state pair 
energy. For small T, the excess density increases as T 2 . 
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(ii) Correction to the leading orde r be havior. For 
T > T*, the right-hand-side of Eq. ( |C12| ) appe ars a s 
a correction to the leading order behavior, Eq. (CIO). 
The divergence of A(T/T*) at T = nT* signals switching 
of behavior for Ii(n), and the corresponding crossover 
can be analyzed in detail by isolating out the contri- 
bution from this term. When the temperature window 
2T* < T < T2 exists, there is another correction term to 



( |CT0| ) from Ix(2) ~ y 2 p {z 2 ) ~ Ny 2 exp(2A/T 2 ). All terms 



The discussion of Appendix C indicates that a possi- 
ble source of singularity comes from the terms I\{n) at 
T — T n , where the argument of the error function in Eq. 
dC9|) undergoes relatively rapid change. A true singular- 
ity, however, appears only in the limit R — > 00. Since 
the statistical weight of large pairs vanish rapidly with 
R, the rapid change of Ii(n) at large R may not produce 
a singularity in F v . The following calculation supports 
this idea. 



h (») = TTtfJ exp(-0„Z) [1 - ert(8 n lV*)] dl (D3) 



where 



other than these two are shown to be of order (W) 2 or In terms of I = ]n(R/a), Eq. (jCSj) can be expressed as 
smaller. Since our treatment of the pair-pair interactions 
is not accurate enough to produce the coefficient of the 
(W) 2 term, these high order corrections will not be con- 
sidered. 

The above analysis shows explicitly that a perturbative 
calculation of (W) in y p is dangerous at all temperatures. 
Even for T > T» , one encounters difficulties when the cal- 
culation is carried out to sufficiently high orders, though 
low order terms are well behaved. Such behavior is typ- 
ical for a function which has an essential singularity at 
V P = 0. 

The calculation of (W 2 ) can be reduced to the calcu- 
lation of (W) with the help of the identity 



= 2-nnK - 4 - 2-nan 2 K 2 , 



S n = 2nT g /T - 2T g /T, 



(D4) 



(D5) 



Substituting (D3) into (Dl), and integrating over /, we 
obtain the contribution to F v from h(n), 



(W 2 



-y, 



dy p 



(C16) 



The right-hand-side of the above equation can be evalu- 
ated using Eqs. (CIO) and (C12) in respective re gimes 
For T < T„, the leading order result is given by ( 4.10 ) 
For T > T«, the leading order expression of (W) is pro- 
portional to y p and hence does not contribute to (W 2 ) 
Going back to Eq. (CI), and keeping the lowest nonvan- 
ishing terms, we obtain, 



(W 2 )^h(2) + B(^) 



NT, 



(27TA) 1 /: 



■ exp 



-—) 

2T?J' 



(C17) 



(D6) 



In deriving the above equation, we have neglected a weak 
dependence of the parameters a and K = J/T on R, 
which is justified asymptotically in the ordered phase. It 



is clear from (D6) that the free energy has no singularity 
at 8 n = or any other point in the ordered phase for any 
n . 



where B(u) = (1 — u)A(u) + u/(u — 2). It is seen that, 
for T > T», (W 2 ) is smaller than (W) by a fact or which 
decreases as a power-law of R. Equation (C17) also ac- 
counts fo r the crossover regime around T = T„, and re- 
duces to d4~l~0|) for T < T*. 



APPENDIX D: VORTEX FREE ENERGY IN THE 
ORDERED PHASE 

Following the same idea as in the calculation of (W) 
in Appendix C, we may rewrite ( 5.22| ) as, 



<IK - -^(h + Y, — [h(n)+h(n)]), (Dl) 



where 



R 2 



I 3 = I (lny p + \nz)P(z)dz. 

'l/Vp 



(D2) 
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FIG. 1. Previously proposed phase diagrams of the disor- 
dered XY model, (a) Order-disorder transition at T = T+ 
and then again at a reentrance temperature T = T_ . (b) No 
reentrance transition, but freezing of vortex-pair excitations 
below T* (dashed line). 
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FIG. 2. Phase diagram of a single vortex. A true glass 
transition takes place at T 9 (solid line). The free energy of 
the vortex vanishes along the dashed line. 
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O jt/4 jt/2 k~ 1 =T/J 

FIG. 5. Schematic phase diagram from present work. 
Solid line indicates the order-disorder phase boundary in 
terms of bare parameters of the model. It lies inside the re- 
gion enclosed by the dashed line, which is the line of phase 
transition when renormalized values are plotted. At the spe- 
cial point S and its counterpart 5", where the T*-line meets 
the phase boundary, a slightly different critical behavior is 
expected. See text. 



FIG. 3. Cayley tree representation of correlations in the 
disorder potential on a single vortex, (a) Division of the dis- 
order potential into subsums over rings centered at vortex po- 
sition i. (b) Representation of the potential by the energy of a 
path on the Cayley tree. The energy of a path is the sum over 
the energies assigned to the nodes, (c) Two sites i and j pick 
up nearly identical contribution from distant quenched dis- 
order, but completely different contribution from inner shells 
(shadowed area) surrounding each site. 





FIG. 6. The RG flow trajectory on the Y 2 -K^ plane. 
(This example corresponds to the lowest of the four flow lines 
illustrated in Fig. 4.) At the transition, the bottom of the 
curve touches the horizontal axis (D — 0). Generically, the 
curve is quadratic around the minimum at K' 1 = KJ 1 . How- 
ever, for KJ 1 — K^ 1 , a cubic singularity is found instead. 
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FIG. 4. Renormalization group flows on the a-K^ 1 plane. 
Solid lines are trajectories of the RG flow, with arrows indi- 
cating flow direction. The flow follows a parabola up to T* 
(dotted line), and then joined smoothly by a horizontal line at 
T > T* . The dashed line is the line of hyperbolic fixed points 
of the RG flow, which gives the phase boundary at Y = 0. 
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